A model for shock wave chaos 
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We propose the following model equation: 

+ ^ (ll^ - UUs),^ = f (X, U^) , 

that predicts chaotic shock waves. It is given on the half-line a; < and the shock is located at a; = 
for any t > 0. Here Us (t) is the shock state and the source term / is assumed to satisfy certain 
integrability constraints as explained in the main text. We demonstrate that this simple equation 
reproduces many of the properties of detonations in gaseous mixtures, which one finds by solving 
the reactive Euler equations: existence of steady traveling-wave solutions and their instability, a 
cascade of period-doubling bifurcations, onset of chaos, and shock formation in the reaction zone. 



Shock waves arise in a wide range of physical phenom- 
ena: gas dynamics, shallow-water flows, supernovae, stel- 
lar winds, traffic flows, quantum fluids, and many others. 
The dynamics of shock waves can be quite intricate and 
difficult to analyze due to the difficult nature of the hy- 
perbolic conservation laws that govern their evolution. 
The theory of shock waves has a rich history beginning 
with the fundamental contributions by Riemann in the 
middle of the 19th century. Nevertheless, numerous open 
questions remain regarding the shock dynamics, espe- 
cially when classical gas-dynamical shock waves interact 
with additional physical or chemical phenomena, such as 
magnetic and gravitational flelds, chemical reactions, ra- 
diation, and others. 

Our focus in this Letter is on one fascinating feature of 
a shock wave propagating in a chemically active medium, 
namely shock-wave chaos. This is a phenomenon wherein 
the shock propagates with its speed oscillating chaoti- 
cally about a certain average. It has previously been 
demonstrated to occur in gaseous detonations, by solv- 
ing the reactive Euler equations [HI [13] . Detonations are 
shock waves in reactive mixtures that are sustained by 
the chemical energy release in the mixture; the reactions, 
in turn, are triggered and sustained by the heating pro- 
vided by the shock compression. 

Analytical and numerical difficulties associated with 
solving the reactive Euler equations motivated the intro- 
duction of simple analog models, with the aim of captur- 
ing the essential nature of observed detonation shocks. 
As it is well-known. Burgers fT] introduced his equation, 
Ut -\- uUx = Uxx (where the subscripts t and x indicate 
partial derivatives) , now a hallmark of hyperbolic differ- 
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ential equations and shock wave theory, in the hope of 
capturing the essential nature of turbulence with a sim- 
ple and tractable model. Following a similar idea, Fickett 
[H |S] and shortly after him Majda |TT] , introduced sim- 
ple analog models for detonations in the hope of gaining 
some insight into the complicated behavior of detonation 
waves, as observed in experiment and numerical simula- 
tions. Fickett's model is a modification of the Burgers' 
equation, which introduces the effects of chemical energy 
release. It takes the form: 

"t + ^("^ + = 0, At=a;(A,u), (1) 

where u is the primary unknown mimicking density, tem- 
perature, or pressure, w is a rate function, and g is a con- 
stant playing the role of a chemical energy release. The 
chemistry here is represented by an irreversible reaction 
reactants products, with A being a normalized con- 
centration of reaction products. At the shock, A = and 
A increases through the reaction zone to reach A = 1 in 
the products. 

Fickett's model has been shown to reproduce some 
of the features of detonations most notably the 

steady-state structures. Still, the key unstable charac- 
ter of detonations had not been reproduced within this 
model until Radulescu and Tang [14] extended it to a 
two-step chemistry with an inert induction zone followed 
by an energy-releasing reaction zone. In [14_ , the authors 
were able to reproduce, with their analog model, the com- 
plexity of chaotic detonations in the Euler equations. 

In this Letter, we propose a model consisting of a single 
equation that predicts steady traveling wave solutions, 
instability through a Hopf bifurcation, and a sequence 
of period-doubling bifurcations with subsequent chaotic 
dynamics. The onset of chaos in our equation appears 
to follow the same scenario as in the logistic map [12J. 
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For the reactive Euler equations and the Pickett's ana- 
log model, the same scenario has been found |H1 [H] ■ We 
note that even though our model is still an analog, it 
is close to the weakly nonlinear model of Rosales and 
Majda [T5], which is rationally derived from the Euler 
equations rather than postulated as the Fickett's (or Ma- 
jda's) models. However, our emphasis here is not on the 
precise relationship of our model to the Euler equations, 
but rather on presenting what we believe is the simplest 
partial differential equation that is capable of capturing 
much of the richness of detonations in the reactive Euler 
equations. 

Our model is the following partial differential equation: 



Ut 



f {X, Us) , 



(2) 



for a; < and t > 0, with an appropriate initial condition, 
u {x, 0). Here Ug (t) = u (0, t) is the boundary value of the 
solution, which is not prescribed but follows by solving 
([2]), as explained below. The source term / needs to 
satisfy certain integrability conditions, as also explained 
further below. 

Equation ([2]) is a simple model for the reaction zone of 
a detonation moving into a uniform state, in coordinates 
attached to the leading shock. This can be seen by ap- 
plication of the Rankine-Hugoniot shock conditions (see, 
e.g. ID) at a; = 0, for ([2| extended by taking / = and 
u = for a; > 0. Indeed, the shock condition 



V[u] + l [u 



1 



[u] = 0, 



(3) 



z+ — z is the jump 
Us and [m^] = —ul), 



where V is the shock speed and [z 
of z across the shock (so that [u] — 
yields ^ = 0. We assume that the shock satisfies the 
usual Lax entropy conditions [5] , so that the characteris- 
tics from both sides of the shock converge on the shock 
That is. 



> and 



dx . 



Us 

2 / x=o- 

Us 



(4) 
(5) 



which require that Us > 0. Therefore, no boundary con- 
dition at a; = is necessary. Finally, we remark that Us is 
a measure of the shock strength, since [u] — ~Us, and for 
that reason we will analyze Us (t) in what follows when 
describing the shock dynamics. 

The most unusual feature of ^ is that the equation 
contains in it the boundary value of the unknown, Us (t). 
This is in fact the key reason for the observed complexity 
of the solutions and has a simple physical interpretation: 
the boundary information from x = is propagated in- 
stantaneously throughout the solution domain, x < 0, 
while there is a finite-speed infiuence propagating from 
the reaction zone back toward the shock along the char- 
acteristics of ([2]). Importantly, in the Euler equations, 
this situation occurs in a weakly nonlinear reactive shock 



wave where the flow behind the shock is nearly sonic rel- 
ative to the shock |15j. One family of acoustic charac- 
teristics is then nearly parallel to the shock, representing 
the slow part of the wave moving toward the shock. The 
second family moves away from the shock and represents 
the influence of the shock on the whole post-shock flow. 
This occurs on a much faster time scale than the infor- 
mation flow toward the shock. Our model makes this fast 
influence instantaneous. 

One can easily obtain the steady-state solution uq (x) 
of ([2]) by solving 



{ul ~ UqUos)' = / ix,uos) 



(6) 



where the prime denotes the derivative with respect to x 
and the subscript s denotes the shock state. The solution 



uo{x)^'^ + J'^ + 2j^f{y,uos)dy. (7) 



The choice of the steady-state shock strength 



uos = 2 




f{y,UQs)dy 



(8) 



corresponds to the Chapman-Jouguet speed in detona- 
tion theory |7], since then the characteristic speed at 
X = —00 is uo (—00) — uos/2 = indicating that the 
sonic point is reached only at an infinite distance from 
the shock. This situation is analogous to the commonly 
used simple-depletion kinetics in gaseous detonations [7j . 
If one substitutes (|8| into Q, the result is 



Uq (x) = 



UOs 




f{y,Uos)dy. 



(9) 



Clearly, for the solution uq (x) to be real and bounded, 
one must require that 



< 



f{y,uos)dy < 00 



(10) 



for any -co < x < 0. This is the constraint on / that 
we mentioned earlier in the Letter. There may be other 
or additional, more stringent conditions on / in other 
circumstances, but their discussion is outside the scope 
of this Letter. 

Now we explore the fully nonlinear and unsteady solu- 
tions of (pi, for the particular case 



1 



2 V4vr/3 



exp 



(a; - Xf jus)) 
4/3 



(11) 



The function / peaks at x/, chosen here as Xf — 
—k (uqs/us)" , where fc > and a > are parameters. 
We first rescale the variables as follows: u by uqs, so that 
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Figure 1: The limit cycles in the plane of the shock strength 
Us (t) vs Us at a = 4.70 (a) and a = 4.85 (b). 



the dimensionless steady-state shock strength is 1. length 
by Z = fc, and time hy t ~ //uos- From Q, putting in 
all the dimensionless variables and rescaling /3 by P, we 
obtain (keeping the same notation for the dimensionless 
variables and parameters) 



Mo (x) = 



1 




1 + crf ((a; + l)/2V;g) 
1 + erf (1/2^) 



(12) 



where erf (x) is the error function. The dimensionless 
form of ^ IS 



4/3 



(13) 



where a = 1/ [4^/4^ (l + erf (l/2V^))] . Equation (13) 
contains only two parameters now, a reflecting the shock- 
state sensitivity of the source function / (an analog of 
the activation energy in detonations) and /3 reflecting the 
width of / (an analog of the ratio between the reaction- 
zone length and the induction-zone length). 

In the computations below, we use the shock-fltting 
algorithm of [8^ on a domain of length L = 10 with 
A'' = 3000 uniformly spaced grid points. We fix /3 = 0.1 in 
all calculations and vary a to capture instability and bi- 
furcations. When long-time data, such as the local max- 
ima of Us (t) are needed, we compute until t = 6000. Sim- 
ulations start with the steady-state solution perturbed by 
numerical noise. In Fig. [T]one can see that a period dou- 
bling occurs as a is increased from ai = 4.70 to a = 4.85. 
Below the critical value ac ~ 4.04, the steady solution is 
found to be stable. 

If a is increased to large values, we observe that 
the reaction zone extends significantly initially, but sub- 
sequently shrinks. Importantly, as the reaction zone 
shrinks, another shock is formed within the reaction zone 
which then overtakes the lead shock at x = 0, exactly 
analogous to what happens in the reactive Euler equa- 
tions. Under these conditions, the dynamics is no longer 
smooth and must be analyzed differently. Therefore, we 
focus on moderately large a, namely a < 5.2 for our 



Figure 2: Values of the local maxima ii™"^ of the shock 
strength as a function of the parameter a; u™"^ are calculated 
at sufficiently large t to make sure the solution has settled to 
its attractor. 




Figure 3: The chaotic attractor in the space of Us, Us, and its 
at a = 5.1. 



particular choice of /3, so that the dynamics is unstable, 
but no internal shock waves appear to form. Remark- 
ably, as a is increased, we observe a sequence of period- 
doubling bifurcations that leads to chaotic solutions at 
a close to or slightly larger than 5, as seen in Fig. [2] 
The onset of chaos apparently follows the same scenario 
as in the logistic map [T^ I17j . The bifurcation diagram 
in Fig. [2] was computed by solving (13 1 until t = 6000 
for the range of a from 3.9 to 5.2, with an increment of 
0.005. For each a, we find the maxima of Ug (t) between 
t = 5000 and t = 6000 and plot them on the figure. 
Based on a sequence of three period doublings, we es- 
timated the Feigenbaum constant ^ |3j to be about 4.5. 
This is in rough agreement with the well-known value of 
6 = 4.669... for the logistic map as well as that found for 
detonations [51 [T? l [Ti l [Tf] . 

We plot the chaotic attractor at a = 5.1 in the space 
of Us, its, and ils as shown in Fig. [3j Its resemblance 
to the Rossler attractor [TB] is evident. Interestingly, 
when we plot the local maxima of Ug versus their prior 
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Figure 4: The Lorenz map showing consecutive local maxima 
(u", of the shock strength Us (t) over large times (from 

t — 3000 to t = 6000) for the chaotic case at a — 5.1. 
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